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ABSTRACT 

We infer the velocity distribution of radio pulsars based on large-scale 0.4 GHz pulsar surveys. We do so 
by modelling evolution of the locations, velocities, spins, and radio luminosities of pulsars; calculating pulsed 
flux according to a beaming model and random orientation angles of spin and beam; applying selection effects 
of pulsar surveys; and comparing model distributions of measurable pulsar properties with survey data using a 
likelihood function. The surveys analyzed have well-defined characteristics and cover ~ 95% of the sky. We 
maximize the likelihood in a 6-dimensional space of observables P, P,DM, \b\,fi,F (period, period derivative, 
dispersion measure, Galactic latitude, proper motion, and flux density). The models we test are described by 12 
parameters that characterize a population’s birth rate, luminosity, shutoff of radio emission, birth locations, and 
birth velocities. We infer that the radio beam luminosity ( i) is comparable to the energy flux of rehtivistic particles 
in models for spin-driven magnetospheres, signifying that radio emission losses reach nearly 100% for the oldest 
pulsars; and (ii) scales approximately as E 1 / 2 * * * * * * which, in magnetosphere models, is proportional to the voltage 
drop available for acceleration of particles. We find that a two-component velocity distribution with characteristic 
velocities of 90 km s -1 and 500 km s _1 is greatly preferred to any one-component distribution; this preference 
is largely immune to variations in other population parameters, such as the luminosity or distance scale, or the 
assumed spin-down law. We explore some consequences of the preferred birth velocity distribution: (0 roughly 
50% of pulsars in the solar neighborhood will escape the Galaxy, while ~ 15% have velocities greater than 1000 
km s" 1 ; (ii) observational bias against high velocity pulsars is relatively unimportant for surveys that reach high 
Galactic \z\ distances, but is severe for spatially bounded surveys; (Hi) an important low-velocity population exists 
that increases the fraction of neutron stars retained by globular clusters and is consistent with the number o 
old objects that accrete from the interstellar medium; (iv) under standard assumptions for supernova remnant 
expansion and pulsar spin-down, ~ 10% of pulsars younger than 20 kyr will appear to lie outside of their host 
remnants. Finally, we comment on the ramifications of our birth velocity distribution for binary survival and the 
population of inspiraling binary neutron stars relevant to some GRB models and potential sources for LIGO. 

Subject headings: pulsars:general — methodsrstatistical 


1. INTRODUCTION 

Only a small fraction of the estimated 10 9 neutron stars (NS) 
in the Galaxy are visible to Earth-bound observers as radio pul- 
sars. In addition to source brightness and survey sky coverage, 
many effects unique to neutron stars determine which objects 
will be detected: (1) highly beamed emission; (2) turning-off 
of radio emission as a NS spins down; (3) migration from a 
thin extreme Population I z distribution at birth to a combined 
thick disk and unbounded population; (4) interstellar plasma 
dispersion and scattering that distort pulses, diminishing search 

sensitivities to distant and rapidly rotating pulsars. Statistical 

studies that aim to determine the birthrate, total number, and 

spatial distribution of radio pulsars must contend with these se- 

lection effects. Indeed, previous studies (some recent examples 

are Lyne et al. 1998; Hartman et al. 1997; Hansen & Phinney 

1997; Lorimer, Bailes & Harrison 1997) differ from one an- 

other significantly in this respect. All analyses have been ham- 


pered by our incomplete knowledge of intrinsic properties such 
as the luminosities and shapes of pulsar beams and their varia- 
tions with spin rate and magnetic field strength. 

We report a new effort to derive the maximum useful infor- 
mation about Galactic neutron stars from the subset observable 
as radio pulsars (excluding millisecond and binary pulsars and 
those in globular clusters), simultaneously solving for many of 
the pulsar properties upon which the effects of survey selec- 
tion depend. We use a new, highly realistic model of NS birth, 
evolution, and detection in radio surveys, and we implement a 
rigorous statistical analysis of the joint distribution of observ- 
able pulsar properties drawn from real surveys. Monte Carlo 
simulations of NS populations are tested for detectability and 
compared, in a 6-dimensional space of observables, to the ob- 
served pulsar sample. We vary the model parameters to infer 
most likely values and credible ranges. To the extent that real- 
world selection effects are adequately described in the model. 
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the analysis yields results that are independent of selection bi- 
ases. 

The main focus of this paper is a determination of the pul- 
sar velocity distribution function at birth and its implications. 
Other results, such as the luminosities of pulsar beams and the 
radio pulsar birthrate, are briefly described here (covariances 
do exist between these model parameters and the best-fit veloc- 
ity results), but a detailed exposition of the model, statistical 
method, and insights into other areas of pulsar phenomenology 
will be presented elsewhere. We summarize the essential fea- 
tures of our model and method in Sec. 2, present our results in 
Sec. 3, and discuss some implications of our new determination 
of the pulsar birth velocity distribution in Sec. 4. 


2. MODEL AND METHOD 
2.1. Summary of the Model 

We model the birth properties and rotational, kinematic, and 
luminosity evolution of a Monte Carlo population of neutron 
stars. The characteristics of eight real-world surveys are then 
used to assess the detectability of each simulated NS as a radio 
pulsar. The model-derived population is compared to the actual 
survey results using a likelihood analysis. 


2.1.1. Neutron Star Birth 


We assume that the NS birthrate is constant in time, uni- 
form per unit area for a Galactic disk with 15 kpc radius, and 
exponentially distributed above and below the disk midplane 
with parameterized scale zo . The orientations of stellar spin 
and magnetic axes are independent and uniformly distributed 
over the sphere. The birth velocity is distributed isotropically 
in direction about the local standard of rest; it incorporates both 
an asymmetric “kick” impulse from a supernova explosion and 
the residual contribution from a disrupted binary, for those ob- 
jects bom in binaries. We do not separately investigate the two 
phenomena here nor do we attempt to model neutron stars that 
remain bound in binaries (high-field pulsars in binaries are a 
negligible population). Two model distributions of the velocity 
magnitude are considered, consisting of one or two Gaussian 
components. The two component model is 
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parameterized by cr V] and cr V2 , the dispersions of the low- and 
high-velocity components, respectively, and by w i , the fraction 
of neutron stars with parent Gaussian of width a Vl . Finally, the 
initial periods and magnetic field strengths are drawn from log- 
normal distributions with parameterized means ((logPo)* for 
Po in seconds, and (logZ? 0 ), for Bo in Gauss), widths (< 7 ] 0gfb , 
^logflo)’ a °d fixed high and low cutoffs. 


2.1.2. Evolution . 

Given initial conditions, the evolution of a new-born pulsar 
to the present epoch is deterministic and described as follows. 
Most of our analysis assumes spin down according to a con- 
stant, magnetic-dipole braking torque 2 with “braking index” 
n - 3. This choice facilitates comparison of our results with 
2 


previous NS population studies, but we have also investigated 
an alternate torque model, with n = 4.5— effects of different as- 
sumptions for n are described in Sec. 3.4. We assume that the 
direction of the magnetic field with respect to the spin axis does 
not evolve with time. Each star moves in the Galactic potential 
given by Paczynski (1990) for the disk and bulge, and by Cald- 
well & Ostriker (1981) for the corona. 

Each star’s luminosity is a function of period and spin-down 
rate, with no intrinsic spread (i.e., pulsars are assumed to be 
standard candles for fixed values of P and P), 

L r = min {p a P^Lo,E} ergs' 1 , (2) 

where (a, /?, Lq) are model parameters, P = P ]5 1CT 15 s s _I , and 
E = 4n 2 IP/P 3 is the spin-down energy loss rate for moment of 
inertia / = 10 45 g cm 2 . Previous analyses have assumed intrinsic 
spreads in the “pseudo-luminosity” (flux x distance-squared) to 
account for viewing of the pulsar beam from different orienta- 
tions as well as possible differences in intrinsic luminosity. We 
choose instead to explicitly model flux variations due to view- 
ing geometry through a beam model, described below, and a 
physically meaningful luminosity (in erg s _1 ). We find that the 
inevitable spread in fluxes due to viewing geometry is substan- 
tial. Some implications of this choice are discussed in Sec. 3.4. 

Slowly rotating neutron stars do not emit in the radio because 
the electric field near the polar cap is too small (e.g., Chen & 
Ruderman 1993). The cutoff may not be abrupt; we implement 
a “death band” as the probability that objects with given P/P 3 
are radio quiet. The probability is 

l/2{tanh [(P l5 /P 3 - lO^)/^] + lj , 

for P in seconds, so that the position and lateral extent of the 
death band are parameterized by DL and o DL . 

2.1.3. Beaming and Pulse profiles. 

For each simulated neutron star, a pulse profile is generated 
using a phenomenological radio-beam model. Our beam model 
was fixed by fitting 0.4 GHz pulse profiles to ten known pulsars 
with viewing geometries that are well-determined from polar- 
ization measurements (e.g., Lyne & Manchester 1988; Rankin 
1 993), and which span much of the P- P region of interest: PSRs 
B0301+19, B0329+54, B0450-I8, B0525+21, B1541+09 
B 1 642-03, B 1 82 1+05, B 1 933+ 1 6, B2002+3 1 , B2 1 1 1+46. Sat- 
isfactory fits were achieved with the following beam descrip- 
tion, which (together with our view that pulsars are standard 
candles for given P , P) is consistent with preliminary results of 
a more ambitious analysis of pulse waveforms and luminosities 
to be presented elsewhere. 

In our model, pulsed radio flux from a NS depends on contri- 
butions from one “core” and one “conai” beam. Each is Gaus- 
sian in form with opening angle p (half-width at 1 je) that scales 
with period as 

Pcore = 1 -S 0 /* -0 ' 5 ( 3 ) 

Pcone = 6.0°P^ 5 . (4) 

The beams subtend solid angles fl = /tt/j 2 / In 2, where the cor- 
rection factor / CODe = 0.4 accounts for the hollowness of the cone 
beam, and /core = 1 * The total radio luminosity /. r is apportioned 
such that the ratio of core to cone peak flux is 

Ecorc / /"cone = P * • (5) 


We omit from our spin-down expression the sin 2 o dependence of the idealized magnetic dipole braking torque, where a is the angle between the 
held axes In magnetosphere models (e.g . Goldreich & Julian 1969), the torque does not vanish for aligned rotators ( n = 0). 
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The core increases in strength with decreasing period and de- 
creasing observing frequency v\ the latter dependence is not 
relevant to the present analysis of 0.4 GHz surveys. An exam- 
ple of the distribution of core and conal flux is shown in Fig. 1. 

Pulse waveforms are constructed for the viewing geometry 
of each Monte Carlo pulsar, with two opjx>sitely-directed radio 
beams centered on the star's magnetic axis. The waveform flux 
at pulse phase 0 corresponding to polar angle 6(<p) to the line 
of sight is 

F(4>) = Fcore exp { -6» 2 /Pcore } + ^cone exp { ~{6 -Of/w]}, (6) 

where 0 and w e are the centroid and effective width of the Gaus- 
sian conal beam. Figure 2 shows typical differential and cumu- 
lative distributions of period-averaged flux density that result 
from our beam model and radio luminosity inferred from the 
likelihood analysis (Sec. 3). The highest fluxes correspond to 
nearly aligned viewing geometries, while the tail at low flux 
corresponds to lines of sight grazing the skirts of the radio 
beam. The two peaks in the differential distribution arise from 
the core (narrower beam, higher average flux) and cone (wider 
beam, lower average flux) components. 

2.1.4. Detection. 

The effects of interstellar dispersion and scattering are ap- 
plied to the pulse profile of each potentially detectable NS, 
where the degree of waveform “smearing" depends on the ap- 
plicable survey instrumentation. The DM and SM values used 
are based on the Taylor & Cordes (1993) model, dithered to 
simulate errors in the DM-distance relationship at the 50% 
level. Pulse profiles are Fourier transformed and tested for de- 
tection by each survey with sky coverage in the direction of the 
simulated pulsar. The minimum detectable flux is computed 
for the antenna gain, receiver and sky temperature, bandwidth, 
dwell time, sample interval, and detection threshold (including 
harmonic summing) applicable to each survey. Detection prob- 
abilities were modified by the ratio of successfully searched 
area to the area contained within the survey boundaries on the 
sky. We do not model specific telescope pointings, but assume 
uniform antenna gain across the region surveyed by each tele- 
scope. In reality, antenna gain drops off away from the pointing 
direction, and survey pointings are generically arranged in a 
close grid with adjacent pointings separated by the half-power 
width of one, roughly Gaussian, telescope “beam.’ 1 The as- 
sumption of constant gain results in an enhanced detection rate 
and, hence, a lower inferred birthrate, which we later correct 
by the ratio of volumes for beams with constant and Gaussian 
cross-sections, / ga in = 1 -4. 

In summary, we draw neutron stars independently from a pa- 
rameterized parent population, evolve them according to the 
simulation model and test for detectability as pulsars in the 
given surveys. This “forward evolution" from birth to detec- 
tion yields distributions of simulated pulsar properties that may 
be compared to those of the observed sample. 

2.2. Surveys and Data : Selection Criteria 

The surveys modeled, all at or near 0.4 GHz, are those 
of Camilo, Nice & Taylor (1996), Dewey et al. (1985), Fos- 
ter et al. (1995), Manchester et al. (1978), Manchester et al. 
(1996), Nice, Fruchter & Taylor (1995), Stokes et al. (1985), 
and Thorsett et al. (1993). Together, they cover some 95% of 
the sky. The data are culled from the May 1995 electronic ver- 
sion of the Taylor et al. (1993; 1995) catalog of 721 pulsars, 


supplemented with 66 P measurements from D’Amico et al. 
(1998), including the 8.5 sec pulsar of Young et al. (1999). Of 
the full set of known pulsars we keep only those that were de- 
tected in at least one of these surveys; for example, we exclude 
objects found in deep searches that targeted supernova remnants 
or X-ray sources, and pulsars detected only in high-frequency 
(1.4 GHz or higher) searches. In addition, to avoid contamina- 
tion from objects with complicated evolutionary histories, we 
exclude from our analysis: binary pulsars, pulsars in globular 
clusters, extragalactic (LMC and SMC) pulsars, and all pulsars 
with spin-down rates less than 10 -16 s s _1 because of potential 
confusion with spun-up objects. The resulting dataset contains 
jVp S r = 435 pulsars, 79 of which have measured proper motions 
or useful upper bounds. An additional 54 cataloged pulsars 
meet the criteria above but do not have measured spin-down 
rates. We account for these objects in our derived birthrate 
through a correction factor f QO p = (435+ 54)/435, but otherwise 
neglect them in the likelihood analysis. 

2.3. Likelihood Function 

The distributions of measured quantities that characterize 
pulsars arise from variations intrinsic to the NS population, bias 
in survey selection, and measurement errors. We account for 
measurement errors as reported in the literature, and model in- 
trinsic variations through Monte Carlo (MC) realizations sub- 
ject to the same detection criteria used to obtain the observed 
pulsar sample. Measurable quantities (“observables") are of 
two types: (a) V a = {/>,/>, |6|,DM}, which are known with suf- 
ficient precision that their distributions are determined by the 
biased sampling of the population; and (b) V b = proper 

motion and flux density, which carry significant measurement 
uncertainty described by observer-reported or assumed likeli- 
hood functions £^ S F} = C^Cf\ We assume a Gaussian form 

for and log-normal distributions above a detection thresh- 
old for Cf s . 

We do not know a priori the form of the likelihood func- 
tion for V a y and therefore construct it numerically. We first 
partition the 4-dimensional space for V a using the real sam- 
ple of pulsars to define bins in a heirarchical fashion — each 
bin corresponds to a small volume element centered on spe- 
cific values of P,P,\b\, and DM. We then create MC pulsars 
according to a specific population/evolution/luminosity model, 
and find the total weight W, of detected objects with observ- 
ables that fall within the boundaries of the i-th bin. For birthrate 
j V, the expected number of such objects is {«,) = N(T M c/^Wi- 
Here, W is the total weight of MC neutron stars bom uniformly 
distributed over time interval 7mc- The likelihood for the V a 
dataset is then formed using products of Poisson probabilities 
for 0 or 1 elements in each bin, 

£,=«^nw. < 7 > 

I 

where the product is restricted to occupied bins, and the total 
expected number of pulsars is 

<"> = £(".)• ( 8 ) 

i 

For precisely measured quantities, the numerical evaluation of 
C„ is exact in the limit of zero bin size and small occupancy 
In practice, we use bins that enclose more than one data pulsar, 
and we generalize the likelihood function accordingly 
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Information from observables must also be included. 
The weight of Monte Carlo points in the i-th bin is W, = 
(where k counts individual MC points); the weight 
of points consistent with the observer-reported uncertainties 
is Ek^ F}k ( sum over joint likelihood of proper motion 
and flux for data in the i-th bin), which may be rewritten as 

w > (j2k Wk ^,F } l / lv >) = w i((£? } ' > ))- The expected num- 
ber of pulsars in the i-th bin consistent with the flux and proper 
motion observations is N(T M c/W)Wj((Cf s )). The likelihood of 
the whole dataset of precisely and imprecisely measured quan- 
tities is then £ = e~^ Hi (M ({£{**))■ 

Maximizing A = ln£ with respect to N yields the birthrate 



where N d =435 is the number of actual pulsars, W d is the weight 
of detected Monte Carlo points, W is the total weight of Monte 
Carlo points and 7 M c is the timescale of the simulation. In this 
paper we analyze the model parameters using the likelihood 
evaluated at the most likely birthrate Aml- 

A iv=£ ln (^) + £ l "<< i: !*'»' on) 

where a constant term has been omitted. 

Given the likelihood function (eq. 10), we employ the tech- 
niques of Bayesian inference to derive plausible ranges for 
model parameters, and to compare models with different num- 
bers of parameters. (For a review of the essentials of Bayesian 
inference, see Gregory & Loredo 1992.) Bayes's theorem states 
that the posterior probability distribution of a model M given 
a dataset T> is proportional to the probability of the data given 
the model, P(T>\M 1 I) — i.e., the likelihood £ — times the prior 
probability of the model P(M\f) y where / represents assump- 
tions within the model. We assume flat priors over finite ranges 
of the model parameters, so that the peak of the posterior is 
just the maximum of £. We derive an uncertainty estimate (or 
“credible region,” in the parlance of Gregory & Loredo) for 
each model parameter from the marginal probability distribu- 
tion for that parameter, i.e., from the likelihood function inte- 
grated over all other parameters. Similarly, a global likelihood 
for each model may be derived by integrating over all model 
parameters — such global likelihoods offer a rigorous means of 
comparison between models that have different degrees of com- 
plexity as reflected in the number of model parameters. Ratios 
of global likelihoods between models are known as “odds ra- 
tios.” 

2.4. Implementation 

For the T) a observables, we defined a 4-dimensional space 
with 4 partitions in each dimension. With this choice, the 435 
known pulsars were distributed into 164 distinct bins, filling 
64% of the observable phase space. We found that the addition 
of a fifth dimension. Galactic longitude /, to the space of ob- 
servables, or the use of two orthogonal components of proper 
motion (pi and /x/,, say) rather than the magnitude, yielded too 
few MC points in some bins. 

We tested our method by generating simulated pulsar data 
“catalogs” according to a particular luminosity and velocity 
model and attempting to recover the model parameters. The 
mock datasets were constructed to resemble our true dataset, 
both in the number of pulsars and in the number and quality of 


proper-motion measurements. We found that the resulting max- 
imum likelihood parameter values differed from their expected 
values by no more than 20%, with the 99% “credible region” 
of the marginalized posterior probability distribution for each 
parameter (see Sec. 2.3) implying uncertainties of 10% for the 
luminosity parameters and logL 0 , 15% for a, and 30% for 
each of the velocity parameters (cr V] ,<7^, vvi). The model that 
described the mock dataset was reliably contained, in a half- 
dozen runs of the simulation, within the volume of the likeli- 
hood function that contained at most 90% of the probability. 
Further details on mock data recovery tests will be presented 
elsewhere. 

For the likelihood analysis, we generated a Monte Carlo 
dataset of nearly 10^ pulsars for 7 mc = 1 Gyr, drawing the 
points from a large range of parameter values; the large value 
of 7 mc was chosen to readily accomodate the oldest (charac- 
teristic age ~ 300 Myr) cataloged pulsars and a range of brak- 
ing models. We tested each point for detection in any one of 
the eight surveys and used the resultant bin occupancies to de- 
rive the likelihood function, £, as defined above. We tested the 
numerical stability of £ for the best models by increasing the 
number of MC points. 

Computations were carried out on the 136-node IBM SP2 
supercomputer at the Cornell Theory Center; evaluation of a 
single set of model parameters typically required 2.5 minutes 
on 48 nodes running in parallel. Computational limitations pre- 
cluded a full 12-dimensional grid search. Instead, we evaluated 
models on 12 -dimensional “crosses,” iteratively repositioning 
the cross in search of a global maximum. We also searched 
over 3- and 5-dimensional sub-grids of the full parameter space, 
exploring in particular the velocity (cr^cr^, wj) and luminos- 
ity (a, /?,£)) parameters. We have found a well-defined global 
maximum in the posterior probability function for models given 
the data. 

3. RESULTS 
3.1. Birth Properties 

Table 1 displays the model parameter values that maximize 
the posterior probability for models given the data (with flat 
priors) over the indicated ranges. The Galactic birthrate of NSs 
visible as radio pulsars resulting from eq. 9 is 

N = /gain x f no p x Nml = (760 yr) l . (11) 

The first two factors in eq. 11 represent the only parts of the 
analysis that are not explicitly modeled by Monte Carlo, i.e., the 
corrections for cataloged pulsars with unmeasured spin-down 
rates and for variations in antenna gain as described above. 
These should not be confused with the more typical beaming 
corrections. Because our model incorporates a specific descrip- 
tion of the radio beam, no post facto beaming correction is 
needed (or allowed), as it is in previous studies. Nonetheless, 
the value of N under-represents the true Galactic radio pulsar 
birthrate. As shown in Table 1, we impose cutoffs on the birth 
magnetic field (B 0 ) distribution corresponding to the highest 
and lowest field strengths in our sample of cataloged pulsars; 
any unseen pulsars lying outside of this range have not been 
accounted for by our N. The value of N will also depend on 
details of the radial distribution of NS birthsites in the Galac- 
tic disk, which we have fixed in our analysis (Sec. 2.1.1), on 
our description of the pulsar beam (e.g., the assumption of cir- 
cularity; Sec. 2.1.3), and on the assumed spin-down law We 
compare our birthrate with previous estimates in Sec. 4.5. 
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The distribution of spin periods at birth in our model is poorly 
constrained by the data. It is consistent with birth at rapid rota- 
tion rates (as dictated by the presence of the Crab and Vela pul- 
sars in our dataset), but the assumed uni-modal (log-normal) 
form for the distribution cannot test the hypothesis of “injec- 
tion” of pulsars bom with long rotation periods (Narayan & Os- 
triker 1990). The birth field-strength distribution, by contrast, 
is very well constrained, with the mean of the log-normal dis- 
tribution roughly consistent with the (present-day) mean field- 
strength for our dataset (10 12 1 G) as expected, given our as- 
sumption of no torque decay. Cut-offs in the Po and Bo distribu- 
tions at low and high values were fixed for all models, as shown 
in Table 1. Finally, the preferred birth scale-height in our model 
is zo ~ 160 pc, consistent with the results of Cordes & Chemoff 
(1998). 

3.2. Radio Luminosity 

The luminosity law with maximum posterior probability is 

L r = P -1 3 f5 > 5 4 10 29 3 erg s" 1 , (12) 

for emissions in the radio band (t/ > 50 MHz) with power-law 
flux spectra. Comparison of our result with earlier statistical 
analyses is made difficult by the fundamental differences be- 
tween our model and previous descriptions of pulsar luminosi- 
ties. Note, however, that for luminosity proportional to mag- 
netic field strength, the exponents would take on the values 
(a,/?) = (0.5, 0.5), while scaling with the spin-down energy loss 
rate would imply ( a , 3 ) = (—3, 1). We find instead that the lu- 
minosity scaling with spin parameters is similar to that for the 
vacuum potential drop across the magnetic polar cap of a neu- 
tron star (e.g., Goldreich & Julian 1969), L r oc E 1 ' 2 • 

L r is comparable in value to particle energy loss rates E p 
calculated in “polar cap” models for the magnetosphere. For 
example, the energy flux in particles derived by Ruderman 
& Sutherland (1975; their Eq. 29) may be written as E p = 
P -171 /^ 43 10 300 erg s“ l , while Arons & Scharlemann (1979) 
construct a “diode” model which implies E p — P -1 45 Pf > j 4 10 31 2 
erg s' 1 . Comparison of L r with E p suggests that the efficiency of 
radio emission alone (i.e., not including X- and 7-ray emission) 
approaches unity as pulsars age. A more model -independent 
statement can also be made: it is conceivable that the death 
band is related to an upper bound in the efficiency with which 
spin energy is converted to radio emission. If we write L r = e r E 
(e r < 1), our luminosity law implies an efficiency of at least 
e r = 2% for pulsars near the point in the P-P diagram where 
they appear to shut off (P ~ 1 s, P15 ~ 0.1); the largest im- 
plied efficiency among the pulsars in our dataset is e r — 30% 
for J2 144-3933, the most slowly-rotating radio pulsar known 
(P= 8.5 s; Young et al. 1999). The requirement that L r < E im- 
plies a lower bound P, 5 > OJOe^P 2 - 8 , a scaling with period 
not unlike the death band we have assumed in our model. 

We have chosen a death-band slope of 3 (log P 15 = 31ogP + 
DL describes the center of the death-band) in order to constrain 
the efficiency of radio luminosity. Theoretical expectations for 
the slope, however, generally lie between 2 and 3 (e.g., Zhang, 
Harding & Muslimov 2000). To investigate possible depen- 
dence of the form of L r on the assumed death-band slope, we 
have evaluated luminosity models (maximizing the likelihood 
over a grid in (a,0,L o ) with all other parameters fixed) for a 
handful of death-band slopes in this range — we find that the 
maximum-likelihood values of the luminosity parameters re- 
main essentially unchanged. We will parameterize the death- 


band slope in future work. For the_assumed P/P 3 line, the 
maximum likelihood model implies DL = 0.5, and o DL = 1.4. 

Finally, we note that the high-energy (ha > 1 eV) luminosi- 
ties of seven pulsars detected with the EGRET 7-ray telescope 
also seem to scale as E 1/2 (Thompson et al. 1999 and references 
therein). 

3.3. Birth Velocities 

The maximum likelihood one- and two-component birth ve- 
locity distributions are displayed in Fig. 3, together with con- 
tours of log-likelihood evaluated in the vicinity of the peak. 
For the two-component model, ~ 50% of pulsars have 3- 
dimensional velocities greater than 500 km s -1 , and about 15% 
have v > 1000 km s' 1 ; also, 10% have velocities less than 100 
km s~ l . For the single-component model, only a small fraction 
of stars have 3-D velocities greater than 1000 km s _1 or less 
than 100 km s _1 . As is evident in Fig. 3, the two-component 
distribution cleanly differentiates fast from slow pulsars, with a 
deficit of objects at intermediate velocities (300-400 km s’ 1 ) . 
We compare the one- and two-component model likelihoods 
using a Bayesian odds ratio, marginalizing the posterior proba- 
bility over the one and three parameters, respectively, to deter- 
mine global likelihoods for the models. We find that the two- 
component model is favored ~ 10 4 to one. 

Our best one-component velocity model is similar to that 
of Lyne & Lorimer (1994), but the true distribution is appar- 
ently poorly described by a single component. Instead, the 
two-component distribution preferred by our analysis produces 
a significantly larger fraction of high-velocity pulsars (and a 
larger mean velocity, ~ 540 km s" 1 ) than has been required by 
any previous study of the radio pulsar population. We believe 
this is a consequence of ( 1) the unbiased nature of our analysis, 
in contrast to previous efforts that have not accounted fully for 
observational selection effects, and (2) our inclusion of deep, 
high-latitude surveys such as the Parkes 0.4 GHz survey of the 
southern sky (Manchester et al. 1996); see Sec. 3.5 for a de- 
tailed discussion of these effects. 

We have explored models that include a third, very low ve- 
locity (< 50 km s -1 ) component. We find that such models 
produce only marginally improved likelihoods, and acceptable 
models permit at most 5% of NSs in such a low-velocity sub- 
population. In terms of a Bayesian odds ratio, the small im- 
provements in the overall likelihood are insufficient to support 
the added complexity of models that include a third velocity 
component. This result is similar to the finding of Cordes & 
Chemoff (1998) but is more significant because our detection 
model accounts for selection against pulsars deep in the dis- 
persing and scattering medium of the Galactic disk, to which 
low-velocity pulsars would be confined. 

As shown in Table 1, we have also evaluated models that 
assume n = 4.5 and no torque decay. The maximum likeli- 
hood of such models is smaller than the maximum likelihood 
for n = 3 by three orders of magnitude for the same number of 
parameters, implying that canonical dipole spin-down is pre- 
ferred over this one alternative braking model. We note, how- 
ever, that the luminosity parameters ( a,0,L o ) do not change 
significantly, nor does the mean birth magnetic field strength 
(logS 0 [G]). The impact of braking-law assumptions is treated 
in more detail in the following section. 

3.4. Effects of Model Assumptions 
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The maximum likelihood (ML) parameter values shown in 
Table 1 are valid in the context of assumptions made in our 
model — i.e., the quoted uncertainties are purely statistical, re- 
flecting the 68% credible region defined by the likelihood func- 
tion. Two fundamental assumptions may be particularly rele- 
vant to the our birth velocity results: spin-down through mag- 
netic dipole braking, with no field decay, and the standard- 
candle hypothesis of pulsar luminosities. We discuss the po- 
tential effects of these assumptions below. 

3.4. 1 . Dipole (n = 3) braking . 

Because our likelihood function matches multiple simulated 
and observed pulsar properties simultaneously, at least two 
sources of velocity information are available: (1) the distribu- 
tion of present-day Galactic scale heights of radio pulsars to- 
gether with some estimate of their ages, and (2) proper motion 
measurements. We find that the spatial distribution alone sig- 
nificantly constrains birth velocities: if we ignore proper mo- 
tion information, the peak of the posterior probability occurs at 
<r Vl = 80 km s _1 , a V2 = 350 km s -1 , and wi = 0.40, consistent with 
the outcome based on the full likelihood, but with a smaller av- 
erage velocity. Measured proper motions contribute preferen- 
tially to the high-velocity tail of the distribution because some 
middle-aged or old pulsars at low Galactic \z\ can masquerade 
as low-velocity objects until a proper motion measurement re- 
veals that they are moving rapidly — the Guitar Nebula pulsar is 
a relevant example (Cordes, Romani & Lundgren 1993). 

Because spin-down governs the timescale on which neutron 
stars leaving the Galactic disk remain active as radio pulsars, 
the choice of braking index must affect the ML velocity model. 
We may qualitatively assess the effects of different spin-down 
laws and field decay as follows. Braking index measurements 
for the Crab and two other young pulsars yield values n < 3; 
if this holds also for older objects, spin-down proceeds more 
slowly than we have modeled and, for a given pulsar, a lower 
characteristic birth velocity would be needed to reach the same 
Galactic \z\ -height. Field decay, however, works in the oppo- 
site sense: if pulsar luminosities decay with decreasing field 
strength, as seems likely, field decay shortens the time in which 
an active pulsar must travel from its birthplace to its present- 
day vertical distance from the Galactic disk, and larger veloc- 
ities than those we have derived would be necessary. The for- 
mal ML result in the analysis of Cordes & Chemoff (1998) for 
braking index and field-decay timescale tb suggests a ~ 2.5 and 
tb — 6 Myr, or (with similar likelihood) n ~ 4—5 and no field de- 
cay. The values of n and r B in the first case would tend to alter 
our result in competing directions, while the rapid spin-down 
implied in the second case further enhances the high-velocity 
fraction in the derived velocity distribution, as demonstrated by 
our best model with n = 4.5 (the mean velocity increases from 
540 km s -1 to 660 km s -1 ). Furthermore, the proper motion 
measurements that boost the high-velocity component are not 
directly affected by details of the chosen spin-down law. We 
believe, therefore, that our primary conclusion, namely that a 
large fraction of Galactic NS are bom with high velocity, is ro- 
bust. 

3.4.2. Pulsars as standard candles. 

A unique aspect of our simulation is the description of ob- 
served pulsar fluxes as an explicitly modeled consequence of 
viewing orientation across a well-defined radio beam with in- 
trinsic luminosity that depends only on period and spin-down 


rate, with no intrinsic spread. Past practice has been to allow for 
a range of observed fluxes (for fixed distance) from an ad-hoc 
distribution intended to mimic the combined effects of beaming 
orientation and any intrinsic spread in luminosity. Our explicit 
beaming model is a physically well-motivated improvement on 
past practice, and the standard-candle (for given spin parame- 
ters) hypothesis is the simplest allowed by the data. We inves- 
tigate the possible consequences of such a fundamental change 
in the treatment of pulsar luminosities as follows. 

In the top panel of Fig. 2, we show the distribution of ob- 
served fluxes, for a Vela-like young pulsar, that result from our 
beam and best-fit luminosity model (thick line ) and from the 
smeared pseudo-luminosity distribution of Prdszy nski & Przy- 
bycien (1984) as modified by Hartman et al. (1997; see their 
Eqs. 2 and 3) (thin solid line), where a beaming fraction of 3 
has been assumed for the latter, i.e., two-thirds of the simu- 
lated pulsars are assumed to be beaming away from the line of 
sight. Several points of comparison between the two flux distri- 
butions are noteworthy. First, the explicit beam model produces 
observed fluxes spanning a wider range than does the Hartman 
et al. smearing function, a somewhat surprising result that de- 
pends, of course, on the chosen shape of the beam (the outer 
“conal” component falls off as a Gaussian in our model). Sec- 
ond, the proportion of undetectable stars differs between the 
two models by less than a factor of two, a reflection of their 
effective beaming fractions. The difference is primarily due 
to the low-flux tail of the beamed distribution, detectable be- 
tween 10 and 200 mJy — if the beam skirts were to fall off more 
rapidly than a Gaussian, the detection efficiency in our simu- 
lation would decrease, resulting in a higher birthrate through 
Eq. 9 (see Sec. 4.5). Third, when convolved by the Gaussian 
representing measurement error (dotted curve), the shape of 
the beamed distribution begins to resemble that of the distri- 
bution adopted to represent an intrinsic luminosity spread. Fi- 
nally, the two distributions are offset in their absolute flux scales 
by nearly an order of magnitude. Most of this offset is due to 
different scalings of the luminosity with spin parameters: our 
best-fit luminosity law scales as £ 1/2 , whereas Hartman et al. 
assume E ^ 3 . If the latter dependence were adopted for our 
beamed model, the resulting flux distribution would shift to the 
left, approaching the smeared pseudo- luminosity distribution. 

To determine whether our beaming model alone produces a 
flux distribution sufficiently wide to describe the data, we have 
evaluated the likelihood for population models in which the rel- 
evant flux depends not only on the purely geometrical spread 
due to viewing orientation, but also on an additional, intrinsic 
spread in luminosity for objects with a given P and P. We mod- 
eled the intrinsic luminosity by a log-normal distribution with 
mean set by the scale parameter Lo and coefficients a and 0 
as in Eq. 2, and with width equal to e times the mean, where 
we treated 6 as a parameter to be inferred. We examined the 
likelihood function as a function of e and L 0 , holding the other 
parameters fixed. We found that the likelihood decreases mono- 
tonically for all values 0 ^ e ^ 1. Although we have not carried 
out a thorough search of parameter space for a global maximum 
of the likelihood, we interpret this behavior as evidence that ge- 
ometric orientation effects are sufficient, without need for “hid- 
den’ 7 variables, to account for the full spread in observed fluxes 
manifested by the radio pulsar population. 

To provide a degree of separation between our luminosity 
model and the results of the likelihood analysis, we may also 
simply drop flux information from the likelihood function, by 
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ignoring V b observables (Sec. 2.3) and maximizing only the C a 
function, Eq. 7. Assumptions inherent in the luminosity model 
then enter into the simulation only through a given pulsar’s de- 
tectability. We find that the velocity and luminosity parame- 
ter values at maximum likelihood, and the inferred birthrate, 
change little (with the exception, as noted in the previous sec- 
tion, of a decreased average velocity) when flux and proper mo- 
tion are eliminated from the likelihood function. 

We conclude that the impact of our luminosity and beam as- 
sumptions on birth velocity results is likely confined to the ef- 
fects of any insufficiency in the empirically-derived radio beam 
geometry we use, and then only indirectly. Further exploration 
of the effects of different beam properties is warranted and is 
underway. 

3.5. Comparison of Velocity PDF with previous results 

Several attempts to constrain the birth velocities of NSs have 
been described in the literature. (See Cordes & Chemoff 1998 
for a discussion of the analyses of Hansen & Phinney 1997 and 
Lorimer et al. 1997). As emphasized earlier, our method of- 
fers an accounting of selection effects not previously available, 
rigorous treatment of uncertainties in, e.g., proper motion and 
distance, and analysis of the joint distribution of observable pul- 
sar characteristics. Here, we compare our result to those of two 
recent studies. 

Cordes & Chemoff (1998; hereafter CC98) examined the 
kinematics of 49 pulsars younger than 10 Myr with measured 
proper motions. Through a multi-dimensional likelihood anal- 
ysis, but purposely not accounting for selection effects, they 
arrived at a birth velocity distribution for these objects de- 
scribed by <r v , = 175^J km s' 1 , <x V2 = 700^ km s ~‘. and 
wi = 0.86. This distribution, shown in Fig. 3, significantly 
under-represents the high-velocity population implied by our 
analysis. To understand this difference, we have investigated 
the selection effects inherent in the choice of the 49 pulsars an- 
alyzed by CC98: because measurement of proper motion typ- 
ically requires a time baseline of many years, all of the CC98 
pulsars were discovered in some of the earliest pulsar surveys, 
which were substantially less sensitive than present-day surveys 
and directed along the Galactic plane. As a result, the CC98 
pulsars, in addition to being young, constitute a volume-limited 
sample. We find that, if pulsars conform to our best-fit luminos- 
ity law, the detection efficiency of low-latitude, volume-limited 
surveys drops steadily with increasing velocity, so that at 1000 
km s' 1 it is roughly one-third of its low-velocity value. In con- 
trast, recent deep high-latitude (e.g., Foster et al. 1995; Thorsett 
et al. 1993) and all-sky (e.g., Manchester et al. 1996) surveys 
are sensitive to much of the volume occupied by high-velocity 
pulsars (away from the disk where propagation effects limit the 
survey volume). Such surveys combat the velocity-dependent 
selection effects to which early pulsar catalogs were subject. In 
this regard, an indirect observational bias against high-velocity 
pulsars remains unaccounted for in our analysis: by construc- 
tion, our likelihood function assumes that proper motion mea- 
surements are available for an unbiased, representative subset 
of the pulsar catalog. In reality, such measurements have been 
attempted predominantly for nearby, bright pulsars that were 
discovered in early low-latitude surveys. Future proper- motion 
measurements of pulsars discovered in high-latitude surveys 
will mitigate this bias and directly probe the high end of the 
pulsar velocity distribution. 

Fryer, Burrows & Benz (1998; hereafter FBB) simulate 


binary evolution including NS kicks and compare the rela- 
tive numbers of the resulting evolutionary endstates, LMXBs, 
HMXBs, and isolated pulsars, with observation. The analy- 
sis of FBB complements our own because they separate binary 
from kick velocities, but they do not account for selection ef- 
fects associated with pulsar detection. FBB conclude that some 
30% of NS are bom with essentially zero kick, with the remain- 
ing 70% being given large (500-600 kms' 1 ) kicks in addition 
to any orbital velocity retained following binary disruption. Our 
inferred birth velocity distribution is not inconsistent with this 
result: the low- velocity (cr v , ~ 90 km s' 1 ) component could be 
produced by binary break-up. This interpretation leads to other 
difficulties, however, in that the birthrate of binary endstates 
would be much larger than is observed if a significant frac- 
tion of NS (such as our wi = 40%) received zero kick at birth 
(Dewey & Cordes 1987). FBB require only that their simu- 
lation yield birthrates of binaries no smaller than are observed, 
without penalizing models that generate too many binaries. The 
necessity for small kicks in FBB’s analysis is dictated by their 
adopted 1% minimum retention fraction for globular clusters. 
Because most cluster pulsars have different evolutionary histo- 
ries than NSs in the Galactic disk (e.g., as a result of collisional 
interactions and the higher binary fraction within clusters), birth 
velocity distributions relevant to disk pulsars may not apply to 
clusters (see Sec. 4.2). We suggest that the zero-kick fraction 
derived by FBB, 30%, is an upper bound, and that our 90 km s 
component is not due solely to relict orbital velocities. This 
conclusion is supported by the work of De Donder & Vanbev- 
eren (1999), who model the effects of binary membership of 
NS progenitors on the ultimate velocity distribution of single 
NSs, for a variety of binary properties. They find that the distri- 
bution of asymmetry-induced kicks dominates the space veloc- 
ity distribution, with binary effects contributing an unimportant 
excess at low velocities, for any distribution of kicks with an 
average velocity > 150 km s' 1 . If a significant zero-kick com- 
ponent exists (e.g., the distribution suggested by FBB), most 
binaries survive and few low-velocity single pulsars remain. 

4. DISCUSSION 

4.1. Asymmetric Kick Mechanisms 

It is widely believed that the large velocities imparted to neu- 
tron stars at birth arise through a combination of orbital dis- 
ruption and a “kick” reaction to an asymmetric supernova ex- 
plosion. Observational evidence for a kick impulse is found 
in a variety of systems: misalignments between the spin and 
orbital angular momentum axes of the relativistic NS-NS bina- 
ries B 191 3+ 16 and B 1534+ 12 (Wex, Kalogera & Kramer 2000; 
Arzoumanian, Taylor & Wolszczan 1998) and the present- 
day spin-orbit configuration of the binary pulsar J0045-7319 
(Kaspi et al. 19%). Also, direct measurement of the proper 
motions of non-binary radio pulsars yield, in extreme cases, 
velocities far greater than can be provided by orbital motion, 
> 1000 km s' 1 . The Guitar Nebula pulsar, e.g., has a projected 
(two-dimensional) velocity ~ 1600 kms 1 , oriented nearly par- 
allel to the Galactic plane (Cordes et al. 1993). Finally, the sur- 
vival of binary systems into late evolutionary stages, e.g., NS- 
NS, NS-white dwarf, or NS-black hole binaries, depends on the 
magnitude, frequency of occurrence, and preferred orientation 
(if any) of asymmetric kicks. The rates of survival of such sys- 
tems have been estimated, without quantitative knowledge of 
the NS kick velocity distribution, through population synthe- 
sis models including orbital evolution (e g., Dewey & Cordes 
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1987) — such efforts have affirmed the need for a velocity im- 
pulse on timescales shorter than or comparable to the orbital 
motion. 

The kick magnitudes required by observations impose strin- 
gent constraints on theoretical models. No single mechanism 
for producing kicks as large as 1500 kms" 1 has been conclu- 
sively identified; see Lai, Chernoff & Cordes (2001; hereafter 
LCC) for a recent review of proposed mechanisms. Other clues 
to the nature of kicks have been scarce — the timescale on which 
kicks act is not known, and no correlation has been firmly estab- 
lished between the magnitudes or preferred directions of kicks 
and any other pulsar property (Deshpande, Ramachandran & 
Radhakrishnan 1999). A possible exception is suggested by re- 
cent Chandra results for the Vela and Crab pulsars that show 
near alignment of the projected proper motion vectors and spin 
axes for these two young objects (but see Wex et al. 2000 for 
a possible counter-example). LCC explore the implications of 
such a correlation on proposed kick mechanisms. 

For prompt natal kicks, two models have been explored in 
some detail. Large-scale density asymmetries in a presupemova 
core, seeded perhaps by energetic convection in inner burning 
shells driving overstable g-mode oscillations, can produce large 
kicks as a result of non- isotropic shock propagation and break- 
out during the explosion (see Lai 1999 for a review). The role 
of rotational averaging of the momentum impulse, however, 
is unclear — it may be difficult to produce spin-aligned kicks 
in this model (LCC) — and the existence of global asymme- 
tries has not been verified numerically, as computational lim- 
itations have restricted simulations so far to consideration of 
two-dimensional models. Asymmetric emission of neutrinos in 
the presence of a strong magnetic field (due ultimately to CP 
violation in the weak interaction) can also impart momentum 
to a proto-NS, but impulses larger than a few hundred km s' 1 
require very large magnetic fields (10 15 “ 16 Gauss; Arras & Lai 
1999). 

An alternative to large impulsive kicks at birth is the elec- 
tromagnetic rocket effect suggested by Harrison & Tademaru 
(1975), which would accelerate a pulsar with an off-center mag- 
netic dipole gradually over its initial spin-down timescale. This 
mechanism encounters two difficulties. Gravitational radiation 
due to unstable r-mode oscillations (e.g., Anderson 1998) may 
considerably enhance spin-down of a newborn neutron star so 
that the rocket effect acts on a much shorter timescale, resulting 
in a lower final velocity. Also, sufficiently large velocities may 
not be achievable on timescales short compared to the orbital 
periods of the progenitor systems that produce close NS-NS bi- 
naries or systems like PSR J0045-7319. 

The significant fraction of the birth velocity distribution lying 
above 500 km s' 1 cannot have its origin in presupemova orbital 
motion, strong evidence for asymmetric kicks through a mech- 
anism that must be able to produce velocities at least as high as 
1500 km s' 1 . On the other hand, our results also suggest that 
up to 40% of newborn NSs experience only a small, or zero, 
kick. It seems likely, then, that two mechanisms are needed, or 
that the mechanism producing large kicks (e.g., global hydro- 
dynamic asymmetry) acts only in a subset of objects, perhaps 
through a threshold in the explosion energy or the spin rate of 
the progenitor core. As discussed above, kicks are needed to 
avoid producing too-many binary systems containing NSs, sug- 
gesting that some mechanism (e.g., asymmetric neutrino emis- 
sion) responsible for small kicks is active in nearly all instances 
of core collapse. A bimodai birth velocity distribution is, how- 


ever, difficult to understand as a combination of independent 
mechanisms that produce large and small kicks (or a large kick 
and an uncorrelated residual velocity from disruption of an or- 
bit). The convolution of two such independent processes should 
result in a broad uni-modal birth velocity distribution. Our anal- 
ysis demonstrates, however, a clear preference for a distribu- 
tion that can be adequately described (relative to any single- 
Gaussian distribution) as a combination of two Gaussian com- 
ponents with very different dispersions. 

4.2. Escape from the Gravitational Potentials of the Galaxy 
and Globular Clusters 

Leonard & Tremaine (1990) derive a lower bound of V e = 430 
km s' 1 for the escape velocity from the Galactic potential in the 
solar neighborhood; the value of V e is larger at smaller Galactic 
radii. Assuming a typical escape velocity of 500 km s' 1 , we find 
that 50% of pulsars will be bom with sufficient velocity to es- 
cape from the Galaxy for our two-component velocity model, or 
nearly 40% for the one-component model. The larger fraction 
is comparable to that derived by Lyne & Lorimer (1994; here- 
after LL94) based on their one-component distribution. CC98 
derived an unbound fraction of 25%, but argued that this figure 
was likely underestimated by a factor of two due to selection 
against high-velocities in their pulsar sample. 

Drukier (1996) estimates the fraction of newborn NSs re- 
tained by globular clusters as a function of kick velocity, for 
both single stars and binaries containing NSs. It is generally 
believed that clusters must retain at least 10% of the NS bom 
within them, but V e for a globular cluster is just ~ 30 km s' 1 . 
For the LL94 distribution Drukier finds that at most 4% of NSs 
are retained for the most massive clusters, with typical fractions 
much smaller. Comparison of our preferred birth velocity dis- 
tribution with LL94 suggests that the latter significantly under- 
represents the low- velocity population of NS (as also suggested, 
e.g., by Hartman et al. 1997): in our preferred model, the frac- 
tion of objects with birth velocity < 30 km s' 1 is an order of 
magnitude larger than that for the LL94 distribution. To the 
extent that our conclusions for pulsars bom in the disk are ap- 
plicable to NSs in globular clusters, the enhanced low-velocity 
population in our model relative to LL94 should significantly 
increase the cluster retention rate over Drukier’s findings. If 
our low-velocity component includes a significant contribution 
from pre-disruption orbital speeds, i.e., in the limit of zero kick 
for w\ = 40% of the population, the retention fraction is in- 
creased further, perhaps to as large as several tens of percent. 

4.3. Pulsar-Supernova Remnant Associations 

The supernova blast that accompanies NS birth drives a fast 
shock wave in the external medium. Our results for the birth 
velocity distribution of NSs may be used to estimate the frac- 
tion of objects that lie within their host supernova remnants 
(SNR), an important criterion used to assess proposed NS-SNR 
associations. We assume an explosion energy of 10 51 £ 51 ergs 
in a medium of hydrogen density n 0 cm' 3 . For Sedov-Taylor 
expansion (Shull & Draine 1987), the remnant scale at time 
10 4 r 4 yrs is R s = 12.5 (£51 /^o) l ^ 5 t^ 5 pc. A NS with birth ve- 
locity v 3 10 3 km s" 1 passes through the edge of the remnant at 
U = 1 45 (£ 5 ,/rto) 1/3 v 3 V3 . The shock becomes radiative (Shull 
1987) at later times, r sh ~ 1.9 1 El{ H ri? n \0* yr, when the size 
is R sh = 16.2 4V /? pc. Thence, R = R sh (t/t %h ) 2/1 until the 
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remnant slows to about 20 km s 1 and merges with the inter* 

stellar medium. In this regime, t 4 = 1.51 £ 51 n 0 v 3 . 

The CDF for birth velocities less than v is equivalent to the 
CDF for NSs lying within the remnant at time t when v and 
t are related as above for both the Sedov-Taylor and radiative 
regimes. We show these CDFs, modified by the effects of pro- 
jection for a distant observer, in Fig. 4 for several combinations 
of explosion energy and ambient density, as a function of NS 
age for our preferred birth velocity distribution. We find that up 
to 10% of pulsars younger than 20 kyr will appear to lie outside 
of their host remnant shells. Such a small fraction is consistent 
with the set of known and proposed pulsar-SNR associations 
(e.g., Gaensler & Johnston 1995). Of course, to assess the plau- 
sibility of a particular proposed association, age, distance, and 
especially proper motion measurements must be made. 

4.4. LIGO Rates and Binary Inspiral Sources of GRBs 

Both the number and orbital characteristics of binaries that 
survive a supernova explosion depend sensitively on the distri- 
bution of kick magnitudes. A number of studies (e.g., Dewey 
& Cordes 1987; Wijers, van Paradijs & van den Heuvel 1992; 
Brandt & Podsiadlowski 1995; Portegies Zwart & Spreeuw 
1996; Lipunov et al. 1997; Bloom, Sigurdsson & Pols 1999) 
have shown that the birthrates of double neutron star (DNS) 
and low- and high-mass X-ray binary systems decrease with in- 
creasing characteristic kick velocity. For the DNS binaries, the 
resulting orbital period and eccentricity also determine the rate 
at which gravitational radiation is emitted; high-eccentricity 
and short-period systems will typically merge within 10 9 yr, so 
that the rate of merger events will be smaller than the birthrate 
of all DNS binaries. In their final moments, merging systems 
are expected to produce gravitational radiation detectable by 
LIGO. 

Assumed quantities (such as the importance of the high end 
of the initial mass function or the supernova rates in nearby 
host galaxies) contribute some uncertainty to DNS birthrate 
estimates, but the greatest source of uncertainty has been the 
unknown kick velocity distribution, yielding rates spanning in 
some cases more than an order of magnitude (e.g., from 9 
to 384 Myr" 1 per galaxy in the work of Portegeis Zwart & 
Spreeuw 1996). For the highest kick velocities that they model 
(a Maxwellian distribution with 3-dimensional velocity disper- 
sion of 450 km s" 1 , their Model C), Portegeis Zwart & Yungel- 
son (1998) find a DNS birthrate of 17 Myr* 1 and a merger rate 
over 10 Gyr of 12 Myr" 1 . They also find that NS-black hole 
mergers will occur at a rate of 1 per Myr. Similarly, Bloom, et 
al. (1999) find a DNS birthrate of 3 Myr" 1 for their simulations 
incorporating large kicks (drawn from a Maxwellian distribu- 
tion with one-dimensional dispersion of 270 km s" 1 ). The birth 
velocity distribution that we have derived should improve such 
estimates significantly. Belczynski & Bulik (1999) provide an 
expression linking the rate of DNS mergers to the kick veloc- 
ity in one or more Gaussian components based on their binary 
evolution simulation, from which we infer a merger rate of 30 
Myr" 1 . A general feature of two-component kick velocity dis- 
tributions is that the resulting merger rate will be dominated by 
the small-kick component, but systems that remain bound fol- 
lowing larger kicks also contribute through their shorter merger 
times owing to their high orbital eccentricities (see, e.g.. Table 
1 of BPS). The inclusion of a significant fraction of large kicks 
into binary evolution simulations will likely bring their birth 
and merger rate predictions in line with empirical estimates 


based on the small sample of known DNS systems (e.g., Ar- 
zoumanian, Cordes & Wasserman 1999; Kalogeraet al. 2001). 

Models linking the coalescence of compact binaries with (at 
least some) gamma-ray bursts (GRBs) have been developed 
(Pacynski 1986; Goodman 1986), and can account for many 
properties of GRBs, including their energetics, spatial distri- 
bution, and (depending on assumptions about beaming) event 
rate. In this context, the discussion of DNS merger rates above 
applies equally well to GRBs — our inferred birth velocity dis- 
tribution supports published estimates of the merger rate that 
assume the largest kicks. Bloom et al. (1999) model the effects 
of momentum impulses on the systemic motion of a compact 
binary in the gravitational potential of its host galaxy to deter- 
mine the expected range of offsets on the sky between a GRB 
and the host galaxy. They find that the details of the distribu- 
tion of kicks imparted to NSs at birth do not significantly af- 
fect these offsets because large kicks tend to disrupt the binary. 
A system’s survival therefore acts as a velocity filter so that, 
for an ensemble of systems, our larger average kicks will not 
necessarily increase the spatial offsets expected of coalescing 
binaries containing NSs. 

Binary systems containing a stellar-mass black hole are also 
potential sources of gravitational radiation for the LIGO and 
LISA detectors, and may drive certain gamma-ray bursts. The 
systemic velocities of low-mass X-ray binaries containing black 
hole candidates (BHC-LMXBs) appear to be smaller than those 
of NS-LMXBs (White & van Paradijs 1996). Black holes that 
form promptly upon the progenitor star’s collapse are not ex- 
pected to undergo kicks, but some black holes may form after 
passing through a short lived NS (or proto-NS) phase, e.g., if 
accretion of fallback material pushes the compact object’s mass 
over the Chandrasekhar limit (see, e.g., Brandt, Podsiadlowski 
& Sigurdsson 1995; hereafter BPS). Such “delayed formation” 
black holes may receive asymmetric kicks just as neutron stars 
do, with smaller kinematic consequences because of the larger 
mass of black holes. One BHC-LMXB, X-ray Nova Sco 1994, 
does have an unusually high peculiar velocity, ~ 1 10 km s" 1 , 
as do a few high-mass X-ray binaries. As pointed out by BPS, 
formation scenarios of Nova Sco are simplified if one invokes 
an asymmetric kick at birth; moreover, elemental abundances in 
its companion support a supernova origin for the black hole (Is- 
raelian et al. 1999). Extension of our radio-pulsar birth velocity 
distribution to black-hole systems would of course be specula- 
tive, but the properties of Nova Sco suggest that kicks similar 
to those associated with NSs also occur as part of some (per- 
haps rare) black-hole formation mechanisms. Any such kick 
would have implications for the survival rate of binaries con- 
taining a delayed-formation black hole. It should be noted that 
black-hole velocity distributions inferred from objects in binary 
systems may be prone to selection bias against high-velocity 
BHCs, the formation of which may have disrupted the host bi- 
nary. 

4.5. Populations of Isolated Pulsars 

Early expectations that large numbers of old NSs accreting 
from the interstellar medium would be detectable in the soft X- 
ray band have not been borne out by observations (Neuhauser 
& Triimper 1999). The earliest estimates, however, were based 
on an assumed NS velocity distribution peaking at ~ 40 km s" 1 
(Treves & Colpi 1991). The much larger average velocities in- 
ferred by LL94 and CC98 suggested a solution to this discrep- 
ancy, as the accretion rate is expected to scale with velocity 
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as v -3 ; magnetic field decay on timescales less than 10 9 years 
could also explain why so few ISM-accreting NSs are seen (for 
a review, see Treves et al. 2000). As described above, how- 
ever, incomplete accounting for selection effects (those arising 
from radio-wave propagation in the Galactic disk) likely leads 
to an underestimated low-velocity proportion in the analyses of 
LL94, CC98, and others, potentially undermining the conclu- 
sion that the deficit of old accreting NSs is due to higher-than- 
expected characteristic velocities. Our analysis models the ef- 
fects of signal propagation and survey instrumentation, provid- 
ing a more robust census of the low- velocity population. We 
infer a larger low-velocity fraction than do LL94 and CC98, a 
result which, together with our limit on an additional very low 
velocity component (< 5% with birth velocity < 50 km s" 1 ), 
should firmly constrain theoretical expectations for the num- 
ber of nearby isolated NSs accreting from the ISM. As pointed 
out by NeuhSuser & TrUmper (1999), the number of isolated 
NS detections in the ROSAT All-Sky Survey is consistent with 
estimates made by Madau & Blaes ( 1994) based on a low birth- 
velocity distribution modified by diffusive heating of the NSs 
in the Galaxy — the resulting distribution peaks near 90 km s _I , 
similar to the low- velocity component of our preferred birth ve- 
locity model. 

The observational phenomena known as soft-gamma re- 
peaters (SGRs) and anomalous X-ray pulsars (AXPs) are be- 
lieved to be associated with neutron stars. SGRs and AXPs 
are similar in many ways: high (~ 10 14 G “magnetar”) canon- 
ical field strengths, slow rotation (periods in the range 6-12 s), 
and luminosities that are too large to be due to magnetic dipole 
braking alone; alternate sources of energy are field decay and 
accretion of either fallback material from the supernova explo- 
sion or ejecta in the vicinity of the NS (e.g., Chatterjee, Hern- 
quist & Narayan 2000; van Paradijs, Taam & van den Heuvel 
1995). In addition to their distinct high-energy emission char- 
acteristics, a possibly significant observational difference be- 
tween the two types of objects is that most AXPs are found near 
the centers of supernova remnants (so that associations are quite 
certain), while SGRs are found far from SNR centers, mak- 
ing associations less likely and requiring large (~ 1000 kms* 1 ) 
NS velocities if the associations are to be viable (Gaensler et 
al. 2001). Our results suggest that such high velocities are 
not rare among radio pulsars, certainly, and possibly among 
all NSs. Although distinctions between AXPs and SGRs on 
purely kinematic grounds remain controversial, models that in- 
voke episodic interaction (e.g., accretion or spin-down due to 
the propeller effect) with fall-back disks or clumpy ejecta will 
clearly be influenced by the assumed velocities of the young 
neutron stars. 

SGRs and AXPs may represent as much as half of the Galac- 
tic NS population (Gotthelf & Vasisht 1999), yet they are dis- 
tinct from radio pulsars in their spin characteristics and ener- 
getics, with no apparent overlap between the two populations. 
Estimates of the NS birthrate based on radio pulsars, including 
our own, must therefore be revised upward to account for these 
objects. 

Our derived birthrate applies to isolated, rotation-powered, 
radio emitting pulsars, although the radio beam need not be 
visible to Earth-bound observers. Requiring non-binarity for 
our data sample results in the exclusion of just four pulsars 
that otherwise meet all of our selection criteria, suggesting that 
our birthrate essentially accounts also for radio pulsars in long- 
lived, non-interacting binary systems (selection effects that re- 


duce survey sensitivities for binaries are unimportant for long- 
period pulsars). From analyses of the local pulsar population, 
Lyne et al. (1998) derive a Galactic pulsar birthrate of one ev- 
ery 60 to 330 yrs, while Hartman et al. (1997) find model- 
dependent rates that cluster around (300 yr)~ ! . Our birthrate 
estimate is a factor 2.5 lower still. Much of the discrepancy be- 
tween our value and that of Hartman et al. ( 1997), however, can 
be attributed to the different effective beaming fractions in our 
analyses, evident in Fig. 2 (Sec. 3.4), suggesting perhaps that 
the low end of the Lyne et al. (1998) range is more appropri- 
ate. As emphasized earlier, our inferred birthrate depends on 
the assumed beam model and spin-down law (with no field de- 
cay), and cutoffs in the distribution of magnetic field strengths 
at birth. Disagreement between birthrate determinations is also 
likely due to different assumed radial distributions of birthsites 
in the Galaxy (Sec. 2.1.1). It is noteworthy that the major- 
ity of young cataloged radio pulsars were discovered in high- 
frequency (i/ ~ 1.4 GHz) surveys of the Galactic Plane, and 
typically not detected, because of propagation effects, by the 
0.4 GHz surveys we have simulated. A future extension of our 
simulation to include high-frequency surveys will more directly 
probe the population of young pulsars near their birthsites, pro- 
viding valuable constraints on their spatial distribution. 

A number of radio-quiet but apparently rotation-driven (i.e., 
distinct from SGRs and AXPs) young neutron stars have re- 
cently been discovered at X- and 7-ray wavelengths, often as- 
sociated with supernova remnants (e.g., see Brazier & Johnston 
1999 for a compilation); these are most naturally described as 
pulsars whose radio beams do not intercept our line of sight, 
and so such objects would also be included in our birthrate 
estimate. However, Brazier & Johnston (1999) derive a NS 
birthrate ~ (91 yr) _l by counting such objects and known ra- 
dio pulsars within a distance of 3.5 kpc and less than 20 kyr 
old, and assuming a Gaussian radial distribution in the Galaxy. 
Again, our radio pulsar birthrate is substantially lower, raising 
the possibility that some young neutron stars may be truly radio 
quiet, i.e., not emitting radio beams. Such a conclusion, which 
cannot be ruled out at present, would be strengthened if a sig- 
nificant discrepancy between the birthrates of radio pulsars and 
young, radio-quiet NSs persists in future analyses of the Galac- 
tic neutron star population. 

5. SUMMARY 

Together with adjustable models of the birth, evolution, and 
beaming characteristics of the Galactic radio pulsar popula- 
tion, we have simulated eight large-area surveys, and accounted 
for observational selection effects, to infer the birth properties 
(velocity, period, magnetic field strength) of pulsars and their 
present-day radio luminosities. Our results suggest that the 
scaling of radio luminosity with spin parameters P and P is 
similar to the scaling of voltage drop across the magnetic po- 
lar cap in pulsar magnetosphere models, and that old pulsars 
convert most of their spin-down energy into radio emissions. 
For birth velocities, we find that two-component distributions 
are preferred over one-component distributions, and we infer 
maximum likelihood velocities for the two (three-dimensional) 
Gaussian components of 90 km s“' and 500 km s 1 , with the 
population roughly equally divided between the two compo- 
nents. This result strongly supports existing evidence that neu- 
tron stars are subject to “kick” impulses at birth, presumably 
through asymmetric supernova explosions, and that the kick 
mechanism must be able to produce velocities of at least 1000 
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km s' 1 . We find that half of all pulsars born near or outside 
the solar circle will likely escape the Galaxy’s gravitational 
pull. The large kicks implied by our birth velocity distribu- 
tion also have important consequences for the survival rates, 
and final orbital configurations, of binaries in which a neutron 
star is formed, as well as for the plausibility of proposed pulsar- 
supernova remnant associations. 

The catalog data that we have used to constrain our mod- 
els provide spin parameters, fluxes, and Galactic locations for 
the known pulsars, but little in the way of direct velocity mea- 
surements through proper motions and accurate distances. Very 
Long Baseline Interferometry techniques now being perfected 
(e.g., Fomalont et al. 1999, Chatterjee et al. 2001) will provide 
high-quality proper-motion and parallax measurements for pul- 
sars within ~ 1 kpc, a substantial improvement to the current 
dataset. Future statistical analyses similar to the one presented 


here therefore promise important gains in our understanding of 
the velocity distribution of neutron stars. 
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PSR B0833-45 (P = 0.089 s, P 15 = 124.84 ss' 1 , d = 0.35 kpc) 
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Fig 2 The differential {top) and cumulative {bottom) distributions of period -averaged flux arising from: {thick line) our beam and best-fit luminosity models, 

eq. 12 given the period, spin-down rate, and distance of the Vela pulsar, for 10 randomly -drawn viewing geometries; {thin line) the spread in -pseudolurrunosjty’’ 
used by Hartman et al (1997) their Eqs. (2) and (3). The dashed line indicates the pulsar’s cataloged flux; the dotted Gaussian curve depicts the likelihood function 
we use to represent measurement error. See text, Sec. 3.4. For our beam description, the hatched histogram represents viewing configurations dial are undetectable 
by the Parkes 0.4 GHz survey (Manchester et al. 1996); the cumulative distribution indicates that a similar all-sky survey with a 10 mjy flux limit will detect 
two-thirds of the Vela-like pulsars at the 350 pc assumed distance. 
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Fig. 3. Differential and cumulative distributions of birth velocities of NSs. The heavy line represents the maximum-likelihood two-component distribution with 
three-dimensional dispersions of the Gaussian components <*., = 90 km s' 1 and a n = 500 km s’ 1 , with w, = 40% of stars drawn from the low-velocity component. 
The thin line represents the best one-component Gaussian model, with dispersion a, = 290 km s' 1 . These results hold for magnetic dipole braking and no field 
decay. The dashed line represents the best-fit velocity distribution determined by Cordes & Chernoff(!998; a, = 175 km s' 1 , <r„ 2 = 700 km s' 1 and tv, = 86%) 
from a sample of young pulsars, without correcting for selection effects in pulsar surveys. Inset: Likelihood contours for the velocity parameters (all others held 
fixed at their maximum-likelihood values). The contour interval is A In £ = -2 throughout. 
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p IG 4 The fraction of neutron stars seen to be outside a characteristic radius of their host remnants (see text. Sec. 4.3). Some objects are physically outside 

but are projected to be interior. Rs is the shell radius and r? is the fraction of the shell radius considered. Different curves are for different explosion energies, 
ambient interstellar densities, and radius fractions. (1) Heavy line: Fraction outside the shell of a SNR with nominal energy expanding into an ISM of unit density. 

= 1 nn = 1 r? = 1 . (2) Fraction outside half the shell radius of a SNR with nominal energy expanding into an ISM of unit density: 4t - 1 ,«o = 1, V - U) 

Fraction outside the shell of a SNR with nominal energy expanding into an ISM with larger density: /$, = l,no = 5,r? = 1 . (4) Fraction outside the shell of a SNR 
with half nominal energy expanding into an ISM of unit density: «, = 0.5, no = l,xj = 1. (5) Fraction outside the shell of a hypernova with x 10 nominal energy 
expanding into an underdense ISM: £ji = 10, no = 0.5,?? = I. (6) Fraction outside three times the shell radius of a SNR with nominal energy expanding into an ISM 
of unit density: £51 = I , «o = l > *7 = 3. 
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Table 1 

Summary of model parameters. 


Parameter 

Peak posterior probability value 2 

Credible 

Range searched 


n = 3 

n = 4.5 

range 

Min. 

Max. 

a 

-1.3 

-1.6 

± 0.3 

-2 

0.5 

p 

0.4 

0.5 

± 0.1 

0 

1 

logLo 

29.3 

29.1 

dr 0.1 

27.5 

30.5 

DL 

0.5 

0.8 

± 0.3 

-1.2 

1.2 

<?DL 

1.4 

1.4 

± 0.2 

0 

3.2 

(log Pols]) 

- 2.3 

-1.9 

± 0.2 

-2.6 

-1.6 

^logPo 

>0.2 

0.6 

± 0.2 

0.1 

0.7 

log Pq (Note b) 

-4 

-4 



_ 


log/^ 

0 

0 






(log Bo [G]> 

12.35 

12.35 

± 0.10 

12.00 

12.50 

^log Bq 

0.40 

0.65 

± 0.05 

0.25 

0.50 

logB>° 

11.2 

11.2 

_ 

_ 

_ 

logfl“ 

13.8 

13.8 

_ 

_ 



Zo (pc) 

160 

220 

±40 

80 

250 

Velocity components 0 

One Two 

Two 




cr V| (km s“ l ) 

290 ±30 901$ 

1201$ 


50 

400 

(j\> 2 (km s" 1 ) 

SOO^iso 

9 s 

0 

1 + 

gg 


400 

1400 

Wi 

0.4 ±0.2 

0.45 ±0.2 


0 

1 

In (C/C™,) 

-25.6 0 

-8.1 





a The parameter value for the maximum likelihood (eq. 10) model, and with flat “priors” assumed. 

^ Lower and upper cutoffs of the log-normal distributions for /ft and B 0 were fixed at the values shown. 
c All trial distributions were truncated above 3000 km s' 1 . 



